How numpy handles day-to-day algebra?

import numpy as np

np.__version__

118 False
140 False
umath 8 False
umath 10 False
umath 12 True
core 73 False
Yes!
numeric 29 False
numeric 41 False
numeric 1128 False
numeric 2515 False
numeric 2517 True
core 75 False
core 77 True
core 83 True
core 94 True
142 False
144 True

'1.25.2'

Motivation

In this blog post, we will try to figure out how numpy handles seemingly simple math operations. The motivation behind this exploration is to figure out if there are a few foundational operations behind most of the frequently used functions in numpy. For the sake of right level of abstraction, we will not look into addition, subtraction, multiplication, and division.

Pi

Background

pi is an irrational number and computing its value to a certain precision is a challenging task. This video talks in detail how people used to compute pi in the past. At the time of writing this blog post, Google holds the record for computing pi to the highest precision to 100 trilian digits. They used y-cruncher program (it’s free. try it!) with Chudnovsky algorithm to compute pi. Here are the first 100 digits of pi:

\(3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679\)

Source code

This is how pi is defined in numpy source code upto 36 digits.

#define NPY_PI 3.141592653589793238462643383279502884 /* pi */

Let’s verify it.

print(f"{np.pi:.64f}")

3.1415926535897931159979634685441851615905761718750000000000000000

Hmm, that looks off. From 16th digit onwards, the values are different. Let’s try to figure out why.

pi = 3.141592653589793238462643383279502884
pi = np.array(pi, dtype=np.float64)
pi = f"{pi:.64f}"
np_pi = f"{np.pi:.64f}"
assert np_pi == pi

Okay, so it seems like converting 36 digits of pi to 64 bit precision went wrong from 16th digit onwards. What a waste of last 20 digits of pi due to floating point errors! Anyways, let’s move on.

Power

Let’s find out what happens when you execute the following code in numpy.

number = np.float64(1.1)
number**1.2

1.1211693641406024